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We develop an inhomogeneous mean-field theory for the extended Bose-Hubbard model with a 
quadratic, confining potential. In the absence of this potential, our mean-field theory yields the 
phase diagram of the homogeneous extended Bose-Hubbard model. This phase diagram shows a 
superfluid (SF) phase and lobes of Mott-insulator (MI), density- wave (DW), and supersolid (SS) 
phases in the plane of the chemical potential /j, and on-site repulsion U ; we present phase diagrams for 
representative values of V, the repulsive energy for bosons on nearest-neighbor sites. We demonstrate 
^ ', that, when the confining potential is present, superfluid and density-wave order parameters are 

j^j ■ nonuniform; in particular, we obtain, for a few representative values of parameters, spherical shells 

(3JQf of SF, MI, DW, and SS phases. We explore the implications of our study for experiments on 

cold-atom dipolar condensates in optical lattices in a confining potential 

a 

j3 ■ PACS numbers: 05.30Jp, 67.85.Hj, 73.43Nq 

gi 

I. INTRODUCTION 

Experimental studies of quantum phase transitions [D-Q hi systems of cold atoms in traps, with an imposed optical 
lattice, have led to a renewal of interest in theoretical studies of lattice models of interacting bosons . Examples 
of such transitions include one from a superfluid (SF) to a bosonic Mott-insulator (MI). This transition was predicted 
by mean-field theories, such as those of Refs. [H, 0], and obtained in Monte-Carlo simulations @ of the Bose-Hubbard 
model before it was realized experimentally. 

In addition to the optical-lattice potential, a confining potential, most often quadratic, is present in all cold-atom 
experiments. This leads to inhomogeneities in the phases that are obtained: simulations [1, Q of the Bose-Hubbard 
model, with such a confining potential, and experiments [ToL [Tl| on interacting bosons in optical lattices, with a 
confining potential, have both seen alternating shells of SF and MI regions. 

Mean-field theories for the Bose-Hubbard model were first developed for the homogeneous case [H, 0, [HI • These 
theories were then extended to the inhomogeneous case [lH to develop an understanding of the Bose-glass phase in 
the disordered Bose-Hubbard model. In recent work [l4[ we have shown how the effects of such a confining potential 
can be treated, at the level of mean-field theory, as was done in the Bose-glass case in particular, we have 
provided a natural framework for understanding alternating SF and MI shells, mentioned above. Here we extend 
this inhomogeneous mean-field theory to account for the different types of phases, SF, MI, density- wave (DW), and 
supersolid (SS), which can occur in the extended Bose-Hubbard model fl5j |. 

The principal motivation for undertaking such a study of the extended Bose-Hubbard model is provided by the 
experiments that have obtained a dipolar condensate of 52 Cr atoms [l6|. To understand these experiments we must 
study lattice models of bosons with long-range interactions [TtJ and not merely the Bose-Hubbard model with a 
repulsive interaction between bosons on the same lattice site. The simplest model that goes beyond such onsite 
interactions is the extended Bose-Hubbard model, which allows for repulsive interactions between bosons on nearest- 
neighbor sites and the aforementioned onsite interaction. In addition to SF and MI phases of the Bose-Hubbard 
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model, this extended model can have a density wave (DW) phase, in which the mean density of bosons is different on 
the two sublattices of the hypercubic lattices we consider, and a super-solid (SS) phase (see, e.g., Refs. [TBI. Il8l. [l9|). 

Before we present the details of our study, we summarize our principal results. We first develop a mean- field theory 
for the homogeneous, extended Bose-Hubbard model by developing on the work of our group on Bose-Hubbard models 
for the spinless and spin-1 cases @,[l2[; this yields the SF, MI, DW, and SS phases and the transitions between them, 
which have been studied by a Gutzwiller-type approximation fl5l ] that is akin to, but not the same as, our mean- field 
theory. We then develop an inhomogencous mean-field theory for the inhomogeneous extended, Bose-Hubbard model 
by generalizing our inhomogeneous mean- field theory for the Bose-Hubbard model [l4j . In particular, when we use a 
quadratic confining potential in three dimensions (3D), our theory yields inhomogeneous phases with spherical shells 
of SF, MI, DW, and SS states. The precise way in which these phases alternate depends on the parameters of the 
model; we study a few illustrative cases explicitly for which we present order-parameter profiles and their Fourier 
transforms. We also discuss the experimental implications of our work. 

The remaining part of this paper is organized as follows. In Sec. |TT] we introduce the inhomogeneous extended 
Bose-Hubbard model and then develop an inhomogeneous mean-field theory for it. In Sec. IIIII we present the results 
of our mean-field theory. Section llVI contains concluding remarks; here we give a brief comparison of our work with 
earlier studies and we explore the experimental implications of our study. 



II. MODEL AND MEAN-FIELD THEORY 

We study the inhomogeneous, extended Bose-Hubbard model that is defined by the Hamiltonian 

- = -1 J2 (a\ aj + h.c) + l — ^niifk-l) 

<i,j> i 

+ J t ninj-j t J2»&i, (1) 

<i,j> i 

where t is the amplitude for a boson to hop from site i to its nearest-neighbor site j, z is the nearest- neighbor 
coordination number, < i,j > are nearest-neighbor pairs of sites, h.c. denotes the Hermitian conjugate, a\, cii, and 
hi = a\a,i are, respectively, boson creation, annihilation, and number operators at the site i, the repulsive potential 
between bosons on the same site is U , the chemical potential fi^ controls the number of bosons at the site i, and V is 
the repulsive interaction between bosons on nearest-neighbor sites. In the inhomogeneous case, the chemical potential 
is [ii = /i — Vr-Rf , where is the uniform part of the chemical potential, Vt the strength of the harmonic confining 
potential, i?? = Yl n =i -^n(*)) wnere X n {i), 1 < n < d, are the Cartesian coordinates of the site i and d is dimension 
of the hypercubic lattice (we study d = 3 explicitly); the origin is chosen to be at the center of this lattice. Clearly, 
if we set V = 0, we obtain the inhomogeneous Bose-Hubbard model of Ref. (HI, which we follow in our mean- field 
treatment below. In this paper, we set zt = 1, i.e., we measure all energies in units of zt. 

If t = and Vt = 0, the model (TTJ) exhibits (a) MI phases, which have an integral number of bosons per site, or (b) 
DW phases, in which bosons preferentially occupy one of the sublattices, say A, of the bipartite, hypercubic lattices 
we consider; the MI phases are favored at large values of U whereas the DW phases appear if V is large. A variety 
of DW phases are possible; we denote them by DW M/2; here M is the number of atoms per unit cell and 2 denotes 
that the unit cell is doubled, i.e., the length of its side is 2. For example, when t = and Vt = 0, the DW 1/2 phase 
has 1 boson on sublattice A and none on sublattice B (or vice versa); in DW 3/2 there is 1 boson on sublattice A and 
2 on sublattice B (or vice versa). 

If Vt — but t 0, SF or SS phases can be stabilised because the bosons can hop through the lattice. Nonuniform 
states appear when we allow Vt ^ as we show below via our inhomogeneous mean-field theory. 

We now generalize the intuitively appealing mean- field theory of Ref. Q, which has been developed for the ho- 
mogeneous Bose-Hubbard model and then extended to the inhomogeneous case in Refs. [H, [lj|. Our generalization 
introduces order parameters that are capable of distinguishing between DW, SS, SF, and MI phases. Conventional 
mean-field theories introduce a decoupling scheme that reduces a model with interacting bosons or fermions to an 
effective, noninteracting problem, which can be solved easily because the effective Hamiltonian is quadratic in boson 
or fermion operators. By contrast, the mean-field theories of Refs. @, [HI: f° r * ne case V = 0, decouple the hopping 
term in Eq. (fTJ), which is quadratic in boson operators, to obtain an effective, one-site Hamiltonian, which can be 
diagonalized numerically. To generalize this to the case V > 0, we have to decouple the number operators in the 
extended Bose-Hubbard term in Eq. ([1]). In particular, we decouple the first and third terms of Eq. (fTJ) to obtain an 
effective one-site problem, which neglects quadratic deviations from equilibrium values (denoted by angular brackets) . 
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The two approximations we use are as follows: 

hihj ~ {hi)hj +hi(hj) - (hi)(hj); (2) 

here the superfluid order parameter and the local density for the site i are, respectively, ipi = (a^) and pi = {hi), 
respectively. The approximation @ can now be used to write the Hamiltonian ([1]) as a sum over single-site, mean-field 
Hamiltonians Hf IF as follows: 



H MF ^J2 n 



MF 



— - (<• i\ Pi- fi t i X* \ 

~~zT~ = 2zt ' i_ ' 1 ~ ^ * ^ ' 

1 V 

+ + + -ji^iPi ~ PiPili ( 3 ) 

where the superscript MF stands for mean field, and (pi = j Xii ; Pi = Tl&Pi+S; an d ^ labels the z nearest 
neighbors of the site i. This form of the single-site, mean- field Hamiltonian is suitable for the inhomogeneous case 
with V T > 0. 

For the homogeneous case, we note that the hypercubic lattices we consider are bipartite, i.e., they can be divided 
into two sublattices A and B. Each site on the A (£>) sublattice has z nearest neighbors each one of which belongs 
to the B (A) sublattice. Thus, if Vt = 0, tpi = ipA and pi = pa if i € A and ipi = ipB and pi = ps if i G B, whereas 
4>i = ipB and pi = ps if i £ A and (pi = ipA and pi = pa if i £ £>• If we require chemical potentials that are conjugate 
to pa and ps, respectively, we can introduce pi = pa if i € A and fii = fig if i G B; similarly, we can define creation, 
annihilation, and number operators for each sublattice and hence write the mean-field Hamiltonian §3§ as follows: 

U%£^U%r + n Mr. (4) 



MF 



(a A ip B + o-a^Pb) + ^{iPa^b + V'aV'b) 



+ — {n A pB - PaPb) + -—n A {n A - 1) -n A ] (5) 

i 2 zt zt 



n/MF 



-(a B i' A + a B 4'A) + 7:(ipBip*A + V'sVu) 



zt y " TA 2 

+ — (n B pA - PbPa) + 7^~i n B \ n B - 1) -n B - (6) 

t 2 zt zt 

If Vt > 0, we first obtain the matrix elements of %f IF in the onsite, occupation-number basis truncated 
in practice by choosing a finite value for n max , the total number of bosons per site, for a given initial set of values 
for {ipi,pi}. The smaller the values of U and V and the larger the value of p the larger must be the value of 
nmax] for the values of U, V, and p we consider n max = 6 suffices; we have checked this in representative cases by 
carrying out calculations with n max = 10. We then diagonalize this matrix, which depends not just on ipi and ipi+s, 
but also on pi and Pi+s, to obtain the lowest energy and the corresponding wave function, denoted, respectively, 
by E g (ipi,ip l+S ; pi, p l+ s) and ^ g ({ipi}, {pi}); from these we obtain the new order parameters ipi = {^g{{i>i},{Pi}) I 
a t | ^ 9 {{^i} , {pi})) and pi = {^ g ({ipi}, {pi}) | hi | {pi}))- We then use these new values of ipi and pi as 

inputs to reconstruct Vjf F and repeat the diagonalization procedure until we achieve self consistency of input and 
output values to obtain the equilibrium values tp^ q and p^ q (we suppress the superscript eq hereafter for notational 
convenience). [This self-consistency procedure is equivalent to a minimization of the total energy E g ({ipi}, {pi}) = 
E g (i/ji,ipi + s; pi, Pi+s) with respect to ipi and pi .] Given the form of the confining potential, the self-consistent 
solutions for {ipi,Pi} must have spherical (circular) symmetry in the three-dimensional (two-dimensional) case; we 
use this spherical symmetry in obtaining the self-consistent solutions. If Vt = 0, we only need the four order 
parameters ipA, ipB, PA, and ps so the problem of finding self-consistency solutions is much simpler than it is in the 
inhomogeneous case with Vt > 0. In principle, {ipi} can be complex, but we find, as in earlier calculations fH, Il2l - fl4| . 
that the equilibrium solution is such that {ipi} are real. 
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III. RESULTS 

In this Section we present the results of the inhomogeneous mean-field theory that we have developed in the previous 
Section for the extended Bose-Hubbard model ([T]). We expect that the onsite repulsion between bosons is stronger 
than the repulsive interaction between bosons on nearest-neighbor sites, so we restrict ourselves to V < U. We begin 
with phase diagrams for the homogeneous case with Vp = 0. We then investigate order-parameter profiles in the 
presence of the confining potential, i.e., when Vp > 0. 

In Figs. Q] (a) and (b), we present phase diagrams in the (p, U) plane for the extended Bose-Hubbard model (Q]), with 
V T = and (a) V/U = 0.6 and (b) V/U = 0.4, with SF (gray), SS (red), MI (black), and DW (green) phases; the MI 
phases Mil and MI2 have, respectively, one and two bosons per site; and DW 1/2 and DW 3/2 are, respectively, DW 
phases with one and three bosons per unit cell with side 2; we take the lattice spacing of the underlying hypercubic 
lattice to be 1. The SF phase is favored at small values of U. If we hold p fixed at low values and increase U, the 
system first undergoes a transition to an SS phase and then to the DW 1/2 phase. The lobe of the Mil phase appears 
above the DW 1/2 lobe and the encompassing sliver of the SS phase; the next few DW and MI lobes appear as shown 
in Figs. [T] (a) and (b). Note that the red slivers of the SS phases encompass the DW lobes completely. Furthermore, 
the DW and SS phases grow at the expense of the SF and MI phases as V/U increases, as we expect for the extended 
Bose-Hubbard model ([T|). The phase diagrams of Figs. [T] (a) and (b) are qualitatively similar to those obtained by a 
Gutzwiller approximation in Ref . (l5j . 

We obtain these phase diagrams by monitoring the dependence of the SF and DW order parameters on U, V, and 
p. In Fig. [2] we show representative plots of ip a (red dashed line) and tpb (black full line), on sublattices A and £>, 
respectively, versus p for U = 12, Vt = 0, zt = 1, and (a) V/U = 0.6 and (b) V/U — 0.4. We also show representative 
plots of p a (red dashed line) and pb (black full line), on sublattices A and £>, respectively, versus p for U = 12, zt = 1, 
and (c) V/U = 0.6 and (d) V/U = 0.4. We can distinguish these phases from each other by noting the following: in 
the SF phase ipA = ipB > and pa = Pb', in the SS phases tpA, ipB > 0, ipA ^ ipB, and pa ^ Pb', in the DW phases 
ipA = ipB — but pa ^ Pb', in the MI phases ipA = ipB = and pa — Pb — m, a positive integer. 




4 6 8 10 12 14 4 6 8 10 12 14 

u u 

FIG. 1: (Color online) Phase diagrams in the (p, U) plane for the extended Bose-Hubbard model with Vt — and (a) V/U = 0.6 
and (b) V/U = 0.4 showing SF (gray), SS (red), MI (black), and DW (green) phases. Mil and MI2 denote, respectively, MI 
phases with one and two bosons per site; DW 1/2 and DW 3/2 are, respectively, DW phases with one and three bosons per 
cubic unit cell with side 2; the basic lattice spacing is taken to be 1; and zt = 1. 

We now use the inhomogeneous mean-field theory, which we have developed in the previous Section, to obtain 
alternating spherical shells of MI, SF, DW, and SS phases in the 3D, extended Bose-Hubbard model JT]) with a 
quadratic confining potential. We do this by obtaining the order-parameter profiles {ipi,Pi} and also by obtaining 
in-trap density distributions of bosons at different values of U, V, and p. 

In particular, we use a 3D simple-cubic lattice with 256 3 sites and Vr/(zt) — 0.002; and we study the following two 
representative case: (a) p/(zt) = 30 and V/U = 0.6; and (b) p/{zt) = 19 and V/U = 0.4. With these parameters the 
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FIG. 2: (Color online) Plots of the superfluid order parameters tf) a (red dashed line) and ipb (black full line), on sublattices A 
and B, versus p/(zt) for U /(zt) — 12 , Vr = and (a) V/U = 0.6 and (c) V/U = 0.4. Plots of the density-wave order parameters 
p a (red dashed line) and pb (black full line), on sublattices A and £>, versus p/(zt) for U/(zi) = 12 and (b) V/U = 0.6 and (d) 
V/U = 0.4. 




FIG. 3: (Color online) Plots of (a) pi (red dashed line and points) and (b) ipi (black dashed line and points) versus the position 
X along the line Y = Z = for V T /(zt) = 0.002, U/(zt) = 12 , V/U = 0.6 and p/{zt) = 30; the moduli of the one- dimensional 
Fourier transforms, namely, \p^\ and of the plots in (a) and (b) are plotted, respectively, in (c) and (d) versus the wave 
vector hx/ir- 



total number of bosons Nt — 10 6 , which is comparable to experimental values. Furthermore, this choice of parameters 
leads not only to SF shells and two well-developed MI shells (Mil and MI2) but also to two well-developed DW shells 
(DW 1/2 and DW 3/2) and SS shells. 

Before we study this shell structure let us explore some order-parameter profiles. Plots of the order parameters Pi 
(red dashed line and points) and ipi (black dashed line and points) versus the position X along the line Y = Z = 
are shown in Figs. El (a) and (b), respectively, for V T /(zt) = 0.002, U/(zt) = 12 , V/U = 0.6, and fj,/(zt) = 30. These 
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FIG. 4: (Color online) Plots of (a) pi (red dashed line and points) and (b) ipi (black dashed line and points) versus the position 
X along the line Y = Z = for V T /(zt) = 0.002, U/(zt) = 12 , V/U = 0.4 and p/{zt) = 19; the moduli of the one- dimensional 
Fourier transforms, namely, \p^\ and \iph\i of the plots in (a) and (b) are plotted, respectively, in (c) and (d) versus the wave 
vector kx/n- 
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FIG. 5: (Color online) Plots ofp» versus X, along the line Y = Z = 0, with p./(zt) = 1.2, Vr/(zt) = 0.0005, and (a) U/(zt) = 4.2 
and V/U = 0.6, (b) U/{zt) = 5.2 and V/U = 0.6, (c) U/(zt) = 6.2 and V/U = 0.6, (d) U/(zt) = 4.9 and V/U = 0.4, (e) 
U/(zt) = 5.9 and V/U = 0.4, and (f) U/(zt) = 6.9 and V/U = 0.4. 



plots show that the region near X = is an MI2 phase with pi — 2 and ipi = 0. As we move outwards from here 
(either towards X > or X < 0), we emerge into an SF phase with a nonintcgcr value of pi and ipi > 0; note that pi 
and ipi do not oscillate here as functions of X . At slightly larger values of \X\ the system moves into a very narrow 
SS region in which both pi and ipi are oscillating functions of X. If we increase \X\, this SS phase gives way to a DW 
3/2 phase in which pi oscillates as a function of X but ipi = 0. A further increase in \X\ yields another very narrow 
region of the SS phase; this is then followed by a narrow SF region. As we increase \X\ some more, the Mil phase 
is stabilized; this is followed by a very narrow SF region, which is, in turn, followed by a narrow SS region, and then 
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FIG. 6: (Color online) The moduli of the one-dimensional Fourier transforms, namely, \pk\, of the plots of pi in Figs. EJa), (b), 
(c), (d), (e), and (f) are plotted, respectively, in (a),(b), (c), (d), (e), and (f) here versus the wave vector kx/n. 
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FIG. 7: (Color online) Plots of ipi versus X, along the line Y = Z = 0, with p/(zt) = 1.2, Vr/(zt) = 0.0005, and (a) U/(zt) = 4.2 
and V/U = 0.6, (b) U/(zt) = 5.2 and V/U = 0.6, (c) U/(zt) = 6.2 and V/U = 0.6, (d) U/(zt) = 4.9 and V/U = 0.4, (e) 
U/(zt) = 5.9 and V/U = 0.4, and (f) U/(zt) = 6.9 and V/U = 0.4. 



a DW 1/2 regime. This gives way to a very narrow SS region, as \X\ increases even more; finally we enter a small 
region in which the boson density vanishes. Such profiles of pi and ipi imply the shell structure that we explore below. 

It is also useful to obtain a complementary, Fourier- representation picture of the profiles in Figs. [3] (a) and (b), 
because it might be possible to obtain them in time-of-flight measurements (see, e.g., Eq. (44) in Ref. [l|). Three- 
dimensional transforms of the shell structure can be obtained, but they are not easy to visualize; therefore, we present 
the one-dimensional Fourier transforms of pi(X, Y = 0, Z = 0) and ipi(X,Y = 0, Z = 0) with respect to X. The 
moduli of these transforms, namely, \pk\ and 1^1, of the profiles in Figs. [3] (a) and (b) are plotted, respectively, in 
Figs. [3] (c) and (d) versus the wave vector hx/ft- The principal peaks in these transforms occur at kx = (or 2tt) and 
kx = tt; the former is associated with the uniform MI and SF phases; and the latter arises from DW and SS phases 
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FIG. 8: (Color online) The moduli of the one-dimensional Fourier transforms, namely, \ipk\, of the plots of ipi in Figs.[7][a), (b), 
(c), (d), (e), and (f) are plotted, respectively, in (a),(b), (c), (d), (e), and (f) here versus the wave vector hx/ft. 

in which real-space profiles oscillate as explained above. In an infinite system with no confining potential, these are 
the only peaks; however, as we have seen above, the quadratic confining potential leads to shells of MI, SF, SS, and 
DW phases; this shell structure leads to the subsidiary peaks that appear in Figs. [3J (c) and (d) away from kx = 0, it 
and 2tt. 

Analogs of the order- parameter profiles of Figs. [3J (a) and (b) are given in Figs. 2] (a) and (b) , for pi (red dashed line 
and points) and ipi (black dashed line and points), respectively, versus the position X along the line Y — Z = for 
Vr/{zt) = 0.002, U/(zt) = 12 , V/U = 0.4 and p/{zt) = 19. From these plots we see that, in this case, the sequence 
of phases is SS, SF, Mil, SF, a narrow strip of SS, then DW 1/2, another narrow sliver of SS, and finally a region 
with vanishing boson density. The moduli of the one-dimensional Fourier transforms, namely, \pk\ and \ipk\, of the 
plots in Figs. [4] (a) and (b) arc plotted, respectively, in Figs. [4] (c) and (d) versus the wave vector kx/iT. 

From the profiles in Figs. [3] (a) and (b) and Figs. 2] (a) and (b) it is clear that the precise sequence of MI, SF, SS, 
and DW shells depends on the parameters in the extended Bose- Hubbard model ([1]). We illustrate this for other sets 
of parameter values via representative plots, in Figs. [5] (a)-(f), of the density order parameter pi versus X, along the 
line Y = Z = 0, with p/(zt) = 1.2, V T /(zt) = 0.0005, and (a) U/{zt) = 4.2 and V/U = 0.6, (b) U/(zt) = 5.2 and 
V/U = 0.6, (c) U/(zt) = 6.2 and V/U = 0.6, (d) U/(zt) = 4.9 and V/U = 0.4, (e) U/(zt) = 5.9 and V/U = 0.4, and 
(f) U/(zt) = 6.9 and V/U = 0.4, respectively. The moduli of the one-dimensional Fourier transforms, namely, \pk\, of 
the plots of pi in Figs. 0(a), (b), (c), (d), (e), and (f) are plotted, respectively, in Figs. [5] (a),(b), (c), (d), (e), and (f) 
versus the wave vector kx/ft- The corresponding real-space plots of ipi and the Fourier-space plots of \ipk\ are given, 
respectively, in Figs. [7] and [5] 

These SF, MI, DW, and SS shells appear as annuli [l4| in a two-dimensional (2D) planar section Vz through the 
3D lattice, at a vertical distance Z from the center as shown, for Vz=q, U/(zt) = 12. and Vr/(zt) = 0.002, in Fig. [5] 
(a), with V/U = 0.6 and (i/(zt) = 30, and Fig. M (b), with V/U = 0.4 and (i/(zt) = 19. In the former case, the core 
region near X = Y = Z = 0, has an MI2 phase, whereas, in the latter case, this central region is an SS phase. As we 
move radially outward from the center, shells of other phases appear; the sequence of shells in Fig. |H] (a) is the one 
that results from the order-parameter profiles in Figs. |3j (a) and (b); and the sequence of shells in Fig. [9] (b) follows 
from the profiles in Figs. 0] (a) and (b). 

For any 2D planar section Vz we can calculate integrated, in-trap density profiles such as N m (Z), the number of 
bosons in the p = m MI annuli; similarly, we can calculate N p / q (Z) in the p = p/q DW annuli. [Here m, p, and q 
are intergers; e.g., we study m = 1 or m = 2 and p/q = 1/2 and p/q = 3/2.] We can also calculate the remaining 
number of bosons, e.g., N^Z) = Nt~ N m (Z). For the parameter values of Figs. [3] (a) and (b), illustrative integrated, 
in-trap density profiles are plotted versus Z in Figs. 0(c) and (d), respectively. These in-trap profiles show the total 
number of bosons Nt (light blue full lines), the number of bosons in MI2 and Mil regions, N2 (black line in Fig. [S] 
(c)) and iVi (brown dash-dotted lines), respectively, the numbers of bosons in DW 3/2 (light green line in Fig. |9](c)) 
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and DW 1/2 (dark green line in Fig.|9](d)) regions, the numbers of bosons in SS regions (red dashed line), [Nt — ^2] 
(white full line in Fig. [§] (c)), and [Nt — N±] (blue dashed lines). The outermost gray regions contain no bosons. 
Such integrated, in-trap density profiles have been obtained experimentally in Ref. [10| (see, e.g., their Fig. (3)) for 
cold-atom systems with SF and MI phases; therefore, it should be possible to carry out similar experiments on the 
dipolar systems [l7| that have motivated our study. 




Z Z 

FIG. 9: (Color online) SF (white), MI2 (black), Mil (brown), SS (red), DW 3/2 (light green), and DW 1/2 (dark green) annuli 
in the 2D planar section Tz=o (see text) with U/(zt) = 12 and V T /(zt) = 0.002 and (a) V/U = 0.6 and fJ,/(zt) = 30 and (b) 
V/U = 0.4 and fi/{zi) — 19; for these parameter values, the integrated, in-trap density profiles are plotted versus Z in (c) and 
(d), respectively. These in-trap profiles show the total number of bosons Nt (light blue full lines), the number of bosons in 
MI2 and Mil regions, N2 (black line in (c)) and JVi (brown dash-dotted lines), respectively, the numbers of bosons in DW 3/2 
(light green line in (c)) and DW 1/2 (dark green line in (d)) regions, the numbers of bosons in SS regions (red dashed line), 
[Nt — N2] (white full line in (c)), and [Nt — Ni] (blue dashed lines) in (c) and (d). The outermost gray regions contain no 
bosons. 



IV. CONCLUSIONS 



We have developed an inhomogeneous mean-field theory for the phases and order-parameter profiles of the inho- 
mogeneous, extended Bose-Hubbard model JT]) by generalizing earlier studies for the spinless @, EH and spin-1 [l2[ 
Bose-Hubbard models. In the homogeneous case, our theory leads to SF, Ml, DW, and SS phases and phase diagrams, 
such as those of Figs. [1] (a) and (b); these are qualitatively similar to those obtained by a Gutzwillcr approximation in 
Ref. [l5|. In the inhomogeneous case, i.e., with Vr > in the extended Bose-Hubbard model ([1]), our theory lead to 
rich, order parameter-profiles (see, e.g., Figs. [3J HI [5j and [7]). We have also explored the Fourier-space manifestations 
of these profiles, the structures of the shells of SF, MI, DW, and SS phases and the associated integrated, in-trap 
density profiles for representative parameter values. Such shell structure has been explored for cold-atom systems 
that can be modelled by Bose-Hubbard models [Io|, EH I20H26I but not for the extended Bose-Hubbard model. 

To make a detailed comparison of our results with experiments, the parameters of the Bose-Hubbard model must 

be related to experimental ones [l| as follows: ~ = y ^ L2 £- exp(2^/-p^), where E r is the recoil energy, Vb the strength 

of the lattice potential, a s (= 5.45 nm for 87 Rb) the s-wave scattering coefficient, a = A/2 the optical lattice constant, 
and A = 825 nm the wavelength of the laser used to create the optical lattice; typically < Vb < 22E r . If we use 
this experimental parametrization, we scale all the energies by E r . [In this paper, we set zt = 1, i.e., we measure all 
energies in units of zt.] For the extended Bose-Hubbard case, the relation of our model parameters to parameters in 
dipolar systems [l6l Il7| is not straightforward because of the long-range interactions. However, rough estimates can 
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be made as follows: 

t = J w*(v -v,)^-^/ 2 + V l {v)]w{v -v 3 )d\, (7) 

where i and j are nearest-neighbor sites, w are Wannier functions, and V/(r) = J2 a =x y zYa cos 2 (fc Q a) is the optical- 
lattice potential with wavevector k. Furthermore, 

U = U u = J Mr - r ? )| 2 y in t(r - r')Hr' - r 4 )| 2 d 3 r d 3 r' (8) 



and 



V = U <ijy = J \w{t - ri )| 2 ^ nt (r - r')\w{v' - v 3 )\ 2 d\ d 3 r' , (9) 

with 

Vint = D 2 - + 6(v-r'). (10) 



Here -D is the dipole moment, a s is the s-wave scattering constant, and m is the mass. The s-wave scattering constant 
of chromium is |a( 52 Cr)| = (170 ± 39)a and |a( 50 Cr)| = (40 ± 15)a , where a = 0.053 nm [27[. 

We hope that our work will stimulate experiments designed to explore SF, MI, DW, and SS shells in dipolar- 
condensate systems [H, Qj} . 
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